Studies of thermal conductivity in FPU-like lattices 
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The pioneering computer simulations of the energy relaxation mechanisms performed by Fermi, 
Pasta and Ulam can be considered as the first attempt of understanding energy relaxation and 
thus heat conduction in lattices of nonlinear oscillators. In this paper we describe the most recent 
achievements about the divergence of heat conductivity with the system size in Id and 2d FPU-like 
lattices. The anomalous behavior is particularly evident at low energies, where it is enhanced by 
the quasi-harmonic character of the lattice dynamics. Remakably, anomalies persist also in the 
strongly chaotic region where long-time tails develop in the current autocorrelation function. A 
modal analysis of the Id case is also presented in order to gain further insight about the role played 
by boundary conditions. 
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PACS numbers: 63.10.+a 05.60.-k 44.10. +i 

Almost one century ago Pierre Debjie conjec- 
tured that nonlinearities should be considered as 
a basic ingredient for explaining relaxation and 
transport mechanisms exhibited by real solids. In 
fact, the simplest model a lattice of harmonic 
oscillators yields quite an unphysical scenario: 
any initial excitation does not evolve towards an 
equilibrium state, but frequently returns close 
the initial state, while transport is purely ballis- 
tic since any Fourier component transfers energy 
unaltered through the lattice at the sound veloc- 
ity. In the thirties the introduction of nonlin- 
ear terms in a perturbative quantum mechanical 
description, allowed Rudolph Peierls to obtain 
a successful explanation of the thermodynamics 
of solids at very low temperature. Only twenty 
years later Enrico Fermi, John Pasta and Stanley 
Ulam tackled the problem of studying the relax- 
ation and transport properties of a lattice of non- 
linear classical oscillators. As discussed all over 
this issue, their numerical studies have provided 
inspiration for an astonishingly large amount of 
investigations. In this contribution we aim at sur- 
veying the progresses and the still open problems 
concerning heat conduction in the FPU model. 
In particular, we discuss the divergence of heat 
conductivity with the system size in the Id and 
2d versions of the model. 



I. INTRODUCTION 

At the beginning of the 50's one of the first digital 
computers, MANIAC 1, was available at Los Alamos 
National Laboratories in the US. It had been designed 
by the mathematician J. von Neumann for supporting 
investigations in several research fields, where difficult 
mathematical problems could not be tackled by rigor- 
ous proofs 1] . Very soon Enrico Fermi realized the great 
potential of this revolutionary computational tool for ap- 
proaching also some basic physical questions that had re- 
mained open for decades. In particular, MANIAC 1 ap- 
peared as a suitable instrument for analyzing the many 
aspects of nonlinear problems that could not allow for 
standard perturbative methods. In collaboration with 
the mathematician S. Ulam and the physicist J. Pasta, 
Fermi proposed to integrate by MANIAC 1 the dynami- 
cal equations of the simplest model of a crystal: a chain of 
classical oscillators coupled by nonlinear forces described 
by the Hamiltonian 
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where the integer index i labels the oscillators, whose 
displacements with respect to equilibrium positions and 
momenta are qi and pi, respectively. For the sake of sim- 
plicity Fermi, Pasta and Ulam considered the interatomic 
potentials 
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which have been therafter termed "a" and "/3" models re- 
spectively. The corresponding equations of motion were 
implemented by M. Tsingou into a program containing an 
integration algorithm that MANIAC 1 could efficiently 
compute. 

As discussed elsewhere in this Focus Issue, the FPU 
numerical experiment was intended to test how equilib- 
rium (equipartition) is approached by an isolated set of 
nonlinearly coupled oscillators. On the other hand, the 
connection of this problem to the one of transmission of 
vibrational energy must not have been escaped to Fermi's 
physical intuition. Already in 1914, the dutch physicist 
P. Debye had argued that nonlinearity in the interatomic 
forces is necessary for the finiteness of thermal conduc- 
tivity of insulating crystals, i.e. for the Fourier's law 
J = —kVT to hold. Since, despite the many simplifica- 
tions, this basic ingredient is included in the FPU model, 
one could hope to describe this important physical effect 
in a concrete case. 

Furthermore, the measurement of the time scale for 
approaching the equilibrium state, i.e. the "relaxation 
time" of plane-wave excitations, would have provided 
an indirect determination of thermal conductivity. In- 
deed, the most elementary picture of heat conductivity is 
based on the analogy with kinetic theory of gases where 
k = Cv s £/3, C being the heat capacity, v s the sound 
velocity and I the mean free path. In a lattice, heat car- 
riers are phonons (classically the normal modes), and it 
is necessary to take into account that the latter have dif- 
ferent group velocities, «k = <9u>/<9k, depending on their 
wavenumber k. Accordingly, the above expression for k 
generalizes to 

K OC ^ CkVkTk , (4) 
k 

where we have introduced the relaxation time Tk = ^k/^k 
that can be determined by phcnomcnologically including 
all possible scattering mechanisms (anharmonicity, im- 
purities, boundary effects, electrons etc.) that must be 
determined in some independent way. 

In their numerical experiment Q FPU observed that, 
at variance with the original intuition, the energy ini- 
tially fed in one of the low, i.e. long-wavelength, oscil- 
latory mode did not flow to the higher modes, but was 
exchanged only among a small number of low modes, be- 
fore flowing back almost exactly to the initial state, yield- 
ing a recurrent behavior. Despite nonlinearities were at 
work, neither a tendency towards thermalization, nor a 
mixing rate of the energy could be identified. The dy- 
namics exhibited regular features very close to those of 
an integrable system. Thus, in the spirit of the Debye 
argument, this infinite relaxation time would imply an 
infinite conductivity also in presence of nonlinear forces! 

This indirect consequence of Fermi, Pasta and Ulam's 
work received an intuitively appealing confirmation with 
the discovery of Zabusky and Kruskal's soliton. Gen- 
erally speaking, whenever the equilibrium dynamics of 
a lattice can be decomposed into that of independent 



"modes", the system is expected to behave as an ideal 
conductor. Besides the trivial example of the harmonic 
crystal, this applies also to the broader case of integrable 
nonlinear models characterized by the presence of "math- 
ematical solitons" originating the balance of dispersion 
and nonlinearity. The idea that solitons may play a 
role in heat conduction dates back to Toda and has 
been invoked to explain the anomalous behavior of the 
FPU model as a consequence of ballistic transport due to 
solitons of the modified Korteweg-deVries equation (see 
e.g. 0). Thereby, the existence of stable nonlinear exci- 
tations in integrable systems is expected to lead to bal- 
listic rather than to diffusive transport. As pointed out 
in Ref. 0] , solitons travel freely, no temperature gradient 
can be maintained and the conductivity is thus infinite. 

II. THE NONEQUILIBRIUM FPU MODEL: 
EARLY RESULTS 

The original FPU simulation probed the nonequilib- 
rium transient dynamics under the only effect of internal 
forces. When dealing with systems that can exchange 
energy with external reservoirs, one wishes to introduce 
the effect of external temperature gradients by means 
of nonequilibrium simulations too. In the spirit of lin- 
ear response, the two concepts should be to some extent 
related, since the relaxation of spontaneous fluctuations 
rules the system response. 

In the present context, the natural way to proceed 
consists in putting the system in contact with two heat 
reservoirs operating at different temperatures T + and T_ . 
Several methods have been proposed based on both de- 
terministic and stochastic algorithms 0. Regardless of 
the actual thermostatting scheme, after a transient, an 
off-equilibrium stationary state sets in, with a net heat 
current flowing through the lattice. The thermal con- 
ductivity of the chain k is then estimated as the ratio 
between the time-averaged flux J and the overall tem- 
perature gradient (T + — TJ)jL. Notice that, by this 
latter choice, n amounts to an effective transport coeffi- 
cient including both boundary and bulk scattering mech- 
anisms. The average J can be estimated in several equiv- 
alent ways, depending on the employed thermostatting 
scheme. One possibility is to directly measure the energy 
exchanges with the two baths. A more general definition 
(thermostat-independent) consists in averaging 

J = 2 + 9n)Fn , (5) 

n 

that is a suitable microscopic expression appropriate in 
the context of lattices with nearest-neighbour couplings 
0. Here, F n = —V'(q n +i — q n ) is a shorthand nota- 
tion for the force excerted by the n— th on the n + 1-th 
oscillator, while a is the lattice spacing. 

Historically, this approach has been followed some 
years after the implications of the original FPU exper- 
iment were appreciated by the nonlinear physics commu- 
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nity. The first papers on the FPU model under steady 
nonequilibrium conditions date back to the pioneering 
studies of Payton, Rich and Visscher || and Jackson, 
Pasta and Waters ^(J- I n both cases, the Authors con- 
sidered cubic plus quartic potential terms resulting from 
the expansion of the Lennard-Jones potential. To inves- 
tigate the effect of impurities in the crystal, either a dis- 
ordered binary mixture of masses or random nonlinear 
coupling constants [T(i| were considered. It should be rec- 
ognized a posteriori that those very first computer studies 
attacked the problem from the most difficult side. In fact, 
even before the effect of disorder was fully understood 
in harmonic chains, they studied systems where anhar- 
monicity and disorder are simultaneously present. Never- 
theless, those early works have the merit to have showed 
how the interplay of the two ingredients can lead to un- 
expected results that, in our opinion, are still far from 
being fully understood. Indeed, Ref.[|| revealed that the 
simple perturbative picture in which anharmonicity and 
impurities provide two independent (and thus additive) 
scattering mechanisms does not hold. More precisely, 
the Authors found even cases in which anharmonicity en- 
hances thermal conductivity. A qualitative explanation 
was put forward by claiming that anharmonic coupling 
induces an energy exchange between the localized modes, 
thus leading to an increase of the heat flux. Additional 
questions that have been investigated were the effect of 
disorder on the temperature field and the concentration 
of impurities. Besides the obvious finding that disorder 
reduces the value of heat conductivity (for fixed finite- 
chain length), it was noticed an asymmetric behavior be- 
tween the case of a few heavy atoms randomly added to 
an otherwise homogeneous light-atom chain and its con- 
verse. The smaller values of the conductivity observed in 
the former cases were traced back to the larger number 
of localized modes KM . 



Several studies of lattices under energy fluxes as well as 
attempts of designing easy-to-simulate models followed 
these first studies. The reader is referred to Rcf. 5] for 
a more detailed account. Among others it is worth men- 
tioning the example of the so-called ding-a-ling model [ll| 
that was the first example where the validity of the 
Fourier's law was convincingly shown. Besides this neat 
instance (that actually belongs to a different class of Id 
chains which include the interactions with an external 
substrate |l2j). many (sometimes contradicting) results 
and interpretations appeared in the literature. In ret- 
rospective, these difficulties can be traced back to the 
presence of very strong correlations and slow dynamics 
that was unexpected for such simple models and that 
could not be tackled due to computational limits. In the 
next sections, we illustrate how and why the Fourier's 
law breaks down in the Id FPU model. 



III. THE QUASI-INTEGRABLE LIMIT 

A remarkable property of models like (|TJ is the exis- 
tence of two distinct dynamical regimes for the relaxation 
dynamics. For the FPU-model, the existence of an en- 
ergy threshold was identified numerically by Bocchieri et 
al. [l3j and confirmed by the resonance-overlap criterion 
proposed by Chirikov and coworkers [l4j |. Further nu- 
merical experiments (see Ref. |l5j and reference therein) 
showed that there exists a value of the energy density e c 
below which almost-regular behavior significantly slows 
down the relaxation. This originates from the fact that, 
although primary resonances do not overlap, higher order 
resonances can do, yielding a slower evolution in phase 
space on a time scale that is inversely proportional to a 
power of the energy density |l6l |. Conversely, above e c , 
equipartition is rapidly reached during a typical simula- 
tion run. For this reasons, e c has been termed strong 
stochasticity threshold Although the above analysis 
was mainly focused on the FPU-/? model, it has been 
shown that such a threshold generically exists for several 
lattice models in 1 and 2d |l5| . 

A related finding is the dependence of the maximal 
Lyapunov exponent A on the energy density. An ana- 
lytic estimate of A for the FPU-/3 problem has been also 
provided 

A( e )~( e 1/4 l° ie t eC - (6) 
[ e I for e > e c ; 

This implies that the "correlation" time due to the pres- 
ence of chaotic instabilities (i.e. the inverse of A) may 
become exceedingly large for small enough values of e. 

The existence of long-lasting transients and the fact 
that the Lyapunov time-scale becomes very large both 
suggest that also stationary transport properties below e c 
may be strongly affected by the quasi-integrability of the 
dynamics. To investigate this features, we simulated the 
FPU-/? chain in contact with reservoirs at different tem- 
peratures. For computational convenience, the parame- 
ters m, ijj and [3 have been fixed to unity so that the only 
relevant physical parameter is the energy per particle e. 
With this choice the threshold energy is e c ~ 0.1 [l7j . 
The average temperature T = (T + + T_)/2 has been 
chosen to yield an average internal energy of the chain 
below e c . We have used the Nose-Hoover thermostats 
described in detail in Ref. 5]. In order to fasten the 
convergence towards the stationary state, the initial con- 
ditions have been generated by thermostatting each par- 
ticle to yield a linear temperature profile. This method is 
very efficient, especially for long chains, when bulk ther- 
malization may be significantly slow. Simulations of the 
FPU-/3 model with chains of length up to 65536 sites and 
free boundary conditions exhibit a monotonous increase 
of the finite-size conductivity (see Fig. The growth 
is linear at small lengths and crosses over to a slower in- 
crease. This is best seen as a systematic decrease of the 
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effective exponent 



TABLE I: Fitting parameteters, (Eq.QJ 



Cteff(jV) 



d In k 



din TV' 

from a g ~ 1 to a a c ff ~ 0.5 (see Fig. 0). 



(7) 




FIG. 1: Finite-size conductivity for the FPU-/3 model below 
the strong stochasticity threshold for three different choices 
of the boundary temperatures T± = T ± AT/2; AT = 0.01. 
Nose-Hoover thermostat with response times fixed to 0.5; 
each point results from an average over a single trajectory of 
about 10 6 time units. Dashed lines are the best fits according 
to Eq. ©. 




FIG. 2: The effective exponent a e ff of the finite-size conduc- 
tivity for the FPU-/3 model below the strong stochasticity 
threshold. Points are obtained by evaluating the centered 
differences from data in the previous figure. 
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Transport properties can be analyzed by computing 
the power spectrum S(f) of the heat current J at equi- 
librium. We choose to work in the microcanonical en- 
semble and integrated the equations of motion with peri- 
odic boundary conditions and zero total momentum with 
a fourth-order symplectic algorithm 20}. The spectra 
obtained for three different energies roughly correspond- 
ing to the average temperatures in Fig. are reported 
in Fig. |3| By comparing the results obtained for differ- 
ent chain lengths, it is possible to conclude that finite- 
size corrections are negligible in the considered spectral 
range. From the inset, where the logarithmic derivative 



Ses(f) = 



dlnS 



(9) 



is reported versus the frequency /, one can notice that 
the Lorenzian tail observable at high frequencies starts 
to cross over towards a regime characterized by a weaker 
divergence and presumably controlled by the same mech- 
anisms operating in the strong chaotic regime (see next 
section). Although one can see some evidence of a sec- 
ond plateau only for the largest temperature (T = 0.04 
- see diamonds in the inset of Fig. it is possible to 
investigate how the crossover time t c increases upon de- 
creasing the energy density by suitably rescaling the fre- 
quency axis. As a result, it turns out that t c is roughly 
proportional to e -16 , which is not too far from the di- 
vergence e~ 2 of the Lyapunov time-scale in the weakly 
chaotic regime (see Eq. ©). It is thus reasonable to infer 
that the dynamical mechanisms underlying the slow re- 
laxation and the anomalous conduction are basically the 
same at such low temperatures. 



IV. THE STRONG-CHAOS REGIME 



The data reported in Fig. ^ does not a priori exclude 
the possibility of a slow convergence to a constant value. 
Since we rather expect a power law divergence (see next 
section) we tentatively fitted the conductivity data with 
a form (see also [l9|') 



k{N) N N a w 

The results of the nonlinear fit are shown as solid lines in 
Fig. n and the fitted parameter are reported in Table [I] 
Notice that the exponent a is almost independent of the 
temperature. 



In the previous section we saw how the almost-regular 
features of the dynamics may give rise to a breakdown 
of a macroscopic transport law (the Fourier's law in this 
case). On the basis of the previous discussion, it might 
be surmised that above the strong stochasticity thresh- 
old such anomalies should somehow disappear as the dy- 
namics becomes more "mixing" and the relaxation times 
shorter. It thus came as a surprise [2lJ when the anoma- 
lous transport features of the FPU chain were actually 
discovered to persist also in this regime j2^|. The origin 
of this anomalous behavior is twofold. The first cause is 
the reduced dimensionality. Indeed, strong spatial con- 
straints can significantly alter transport properties: the 
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FIG. 3: Power spectra of the flux J as defined in Eq. © for 
the quartic FPU-/3 model with N = 2048 (solid). Data are 
averaged over about 2500 random initial conditions. To min- 
imize statistical fluctuations, a binning of the data over con- 
tiguous frequency intervals has been performed. The curves 
are for energies e = 0.01, e = 0.02, e = 0.04 (bottom to top) 
and have been vertically shifted for clarity. In the inset, the 
logarithmic derivative is reported versus the frequency, after a 
suitable rescaling of the latter quantity. Circles, squares and 
diamonds refer to e = 0.01, 0.02, and 0.04, respectively. 

response to external forces depends on statistical fluc- 
tuations which, in turn, crucially depend on the system 
dimensionality D. This is very much reminiscent of the 
problem of long-time tails in fluids 23] where, for D < 2, 
transport coefficients may not exist at all. The second 
necessary condition is momentum conservation i.e. the 
existence of long-wavelength (Goldstone) modes that are 
propagating and very weakly damped (rfc diverges for 
small k). This latter condition is necessary for these 
anomalies to occur and means that no external (i.e. sub- 
strate) forces must be present. This is precisely the case 
of the cited ding-a-ling model and of other models 
in the same class (the nonlinear Klein-Gordon chain for 
example |l2j). The only remarkable exception to this is 
the coupled rotor chain where, however, different mech- 
anisms are at work [24| . 

Anomalous behaviour means both a nonintegrable al- 
gebraic decay of equilibrium correlations of the heat cur- 
rent J(t) (the Green-Kubo integrand) at large times 
t — > oo and a divergence of the finite-size conductivity 
k(L) in the L — > oo limit. As a results of a series of 
simulation studies p|, it can be stated that for Id lattice 
models one finds 

k(L) oc L a , (J(t)J(0)) oc t-( 1+ V , (10) 

where a > 0, — 1 < <5 < 0, and ( ) is the equilibrium av- 
erage. For small applied gradients, linear-response the- 
ory allows establishing a connection between the two ex- 
ponents. By assuming that n(L) can be estimated by 
cutting-off the integral in the Green-Kubo formula at the 
"transit time" L/v (v being some propagation velocity of 
excitations), one obtains k oc L~ s i.e. a = —5. 

It is thus natural to argue about the universality of 
the exponent a. On the one hand, two independent 
theoretical approaches, self-consistent mode-coupling ap- 



proximation [25l |2(| |2?J and kinetic theory |28|. yield 
a = 2/5. On the other hand, a renormalization group 
calculation on the stochastic hydrodynamic equations for 
a Id fluid [2£j gives a = 1/3. Validation of one of the the- 
ories is still under debate, and the existence of crossovers 
among different scaling regimes has been observed |26| . 
The available numerical data for a range from 0.25 to 
0.44 0, |3(j. As a word of caution, it must be stressed 
that a numerical estimates are indeed challenging. Even 
in the most favorable case of computationally efficient 
models, as the ID gas of hard-point particles with alter- 
nating masses [U, finite-size corrections to scaling are 
sizeable. As a matter of fact, estimates of a as diverse 
as 0.33 [32, H3 and 0.25 34] for comparable param- 
eter choices have been reported. Although in this lat- 
ter instance the anomaly is possibly due to the lack of 
microscopic chaos j32, it is a generic fact that sensible 
results require reaching the limits of present computing 
resources. For instance, in the case of the FPU chain, 
the best estimate sofar (see below) required simulations 
of up O(10 4 ) particles and O(10 8 ) integration steps (plus 
ensemble averaging) |30j. 

Since we are interested in the long-wavelength and 
small-frequency behavior, it is convenient to consider a 
highly nonlinear model in the hope that the asymptotic 
regime sets in over shorter time and space scales. More- 
over, it is advisable to work with a computationally sim- 
ple expression of the force. A reasonable compromise is 
the "infinite temperature" FPU model [3(| 

V(z) = \z A ■ (11) 

This model has no free parameters: since the potential 
expression is homogeneous, the dynamics is invariant un- 
der coordinate rescaling, so that the energy per particle 
e can be set, without loss of generality, equal to 1. 

The power spectra S(f) of J are reported in Fig. 0] 
The long-time tail l|10|) corresponds to a power-law di- 
vergence f s in the low-/ region. By comparing the re- 
sults obtained for different numbers of particles, one can 
clearly see that finite-size corrections are negligible above 
a size-dependent frequency f c (N). By fitting the data in 
the scaling range [f c {N),f s ], where f s ~ 10~ 3 , we find 
S = —0.39(6). These values are consistent with previ- 
ous, less-accurate, findings for similar models, such as 
the standard FPU [H |H| and the diatomic Toda H2 
chains, thus confirming the expectation that they all be- 
long to the same universality class. 

In order to perform a more stringent test of the scaling 
behavior, we again evaluateid the logarithmic derivative 
(jSJ for different frequencies. Since finite-size effects are 
responsible for the saturation of S(f) when / — * 0, / C (AT) 
can be identified (see Fig.0J as the frequency below which 
5 e g starts growing towards zero. Above f c , the quality of 
our numerical data allows revealing a slow but system- 
atic decrease of d c s upon decreasing /, which approaches 
—0.44, a value that is incompatible not only with the 
renormalization-group prediction of Ref. 29], but also 



6 




FIG. 4: Power spectra of the flux J as defined in Eq. (J^J for 
the quartic FPU model (EJ with N = 2048 (solid) and 1024 
(dashed). Data are averaged over 30,000 random initial con- 
ditions. To minimize statistical fluctuations, a binning of the 
data over contiguous frequency intervals has been performed. 



averaged power spectra of the mode amplitudes 
1 N 

* n— 1 

The spectra display a phonon-like peak at a frequency 
which is in excellent agreement with the one computed 
in Ref. j^l (see Fig. EJ): the estimated sound speed is 
v s = 1.34±0. 04, see above. The peaks' linewidths provide 
a measure of the inverse of the relaxation times; they are 
found to scale as a power of the wavenumber (see Fig.0) 
with an exponenent which is very close to the estimate 
5/3 based on arguments of mode-coupling theory plj . 
Notice that in the standard case one would rather expect 
a quadratic law, corresponding to the usual friction term 
of elasticity theory. 



with the result of mode-coupling [25l |27j and kinetic 
|28| theories. Furthermore, convergence seems not fully 
achieved in the accessible frequency range. 

To check the consistency of equilibrium and nonequi- 
librium simulations, one can assume, following the argu- 
ment exposed below Eq. I|10(l , that the finite-size conduc- 
tivity k(L) (with L = Na) is determined by correlations 
up to time r = L/v s , where v s is the sound velocity. 
This means that the frequency / can be turned into a 
length L — v s / f. It might be argued that the absence of 
a quadratic term in prevents a straightforward defi- 
nition of such a velocity in the T = limit; nevertheless, 
it has been shown 36] that an effective phonon disper- 
sion relation at finite energy density can be evaluated for 
model Ull)). yielding v s = 1.308 at e = 1. Using this 
value, we can ascertain that, at least for N > 1000, there 
is an excellent agreement between the two approaches 
(see again Fig.[5j|. 
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FIG. 6: Structure factors for the quartic FPU model (II IP 
with N = 8192 (solid) for the modes of indexes k = 1,2,4, 8 
(left to right). Microcanonical simulations are performed for 
the energy density e = 1 and averaged over an ensemble of 
about 200 initial conditions. 
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FIG. 5: Quartic FPU model: the effective exponent a c g of 
the finite-size conductivity for T+ = 1.2, T_ = 0.8 (full dots), 
compared with the results (— 5 e ff) of equilibrium simulations. 
The two horizontal lines correspond to the theoretical predic- 
tions, 1/3 and 2/5. 



FIG. 7: The wavenumber dependence of the linewidths for 
the quartic FPU model IIH . The linewidths are obtained 
by fitting the peaks of the spectra for different chain lengths 
with a lorenzian lineshape. The dashed line is a power-law 
fit yielding an exponent 1.677 ± 0.025. 



Another way to detect hydrodynamic anomalies is to 
measure the equilibrium structure factors, namely the 
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V. MODAL ANALYSIS 

Conductivity properties are mostly investigated in real 
space because one is interested in the dependence of the 
heat flux on the system size. Nevertheless, the renor- 
malization group analysis has teached us that a (spa- 
tial) Fourier analysis can be very useful especially to un- 
derstand anomalous scaling behaviours. Additionally, in 
the case of integrable systems, a modal analysis allows 
solving the problem by decomposing it into many inde- 
pendent parts; last but not least the effect of boundary 
conditions can be better appreciated. In spite of such ad- 
vantages, very few numerical studies have been devoted 
to investigating heat conductivity in Fourier space for 
nonlinear chains |M l37j|. Here in the following, we sum- 
marize the current understanding and discuss the open 
problems with the help of numerical simulations to have 
some hints about the expected scenario. 

As proposed in Ref. |37j . it is convenient to start from 
the formal energy balance 

4 + Jk + J'k = 0, (13) 

where j£, JjT , and Jj} 1 denote the energy exchanged (per 
unit time) by the A;th mode with the hot, cold heat bath, 
and with the other modes, respectively. By summing 
over the index k, one finds that 

J + + J' = 0, (14) 

where is the total energy flowing from (to) the hot 
(cold) bath; in fact J^k ~ 0> since the energy of the 
system is, on the average, constant. As a result, the total 
heat flux J can be, e.g., decomposed in modal contribu- 
tions, 

k 

So far, this is just a formal statement. In order to make 
it operative, it is necessary to give an explicit definition 
of the modal fluxes. Using the definition (|12fl . one can 
show that the dynamics of the fcth mode is described by 
the following equation 

Q k + u 2 k Qk + <f>t + <f>t + 4>k = 0, (16) 

where uij- is the frequency of the k mode (which follows 
from an energy dependent renormalization of the har- 
monic frequency), while (j) 1 ^ 1 is the effective force due to 
the interactions with all other modes, while (fy^ account 
for the interaction with the thermal baths. 

Upon multiplying Eq. (|16[) by Qk, and introducing the 
modal energy Ek = (Q\ + u>lQ 2 )/2, it follows that 

= (E k ) = (4>l l Qk) + (cp+Qk) + (^Qk) (17) 

where (•) denotes a time average and the first equality is a 
trivial consequence of stationarity. Since we are assuming 
that the lattice spacing is a = 1, it is now clear that the 



three terms in the r.h.s. of the above equations coincide 
with the modal fluxes JJ} 1 , , and JjT , respectively. 

All such contributions can be easily computed when 
the interaction with the thermal bath amounts to instan- 
taneous stochastic collisions occurring to either the first 
or last particle. In such cases, the fluxes can be eas- 
ily computed by a mode expansion of the energy locally 
exchanged in each collision, while J£ l can be obtained 
from the energy variation in between collisions. 

In a linear chain, the eigenmodes are independent of 
each other, so that the term J£ l vanishes. In such condi- 
tions, the effective equation of the fcth wavevector reduces 
to 

Qk + ujlQk + itQk + Ik Qk + e + £" = 0, (18) 

where 7^ gauges the interaction with the thermal bath 
while £j. is a Gaussian white noise, whose diffusion con- 
stant is, according to the fluctuation-dissipation theorem, 
= fcsT ± 7fe. The above equation is nothing but a 
Langevin equation describing an "oscillator" with a dissi- 
pation 7^+7jT and an effective noise corresponding to the 
intermediate temperature T = (7^7 1+ + ^T^) / '{"ft + 
7jT). Although the average energy of the oscillator nei- 
ther increases nor decreases, the energy exchanged per 
unit time with each bath 

4 = -M) (t + €) (19) 

is different from zero. 

Let us illustrate this under the simplifying but other- 
wise general assumption of equal coupling strength with 
the heat baths (7^ = 7^ = Jk)- In fact, in this case, T = 
T a = (T + +T~)/2 so that the energy exchanged with the 
hot heat bath is J+ = -2j k T + 2 7fc T+ = j k (T+ - T~), 
i.e. Jk ~ 7fe- Conductivity properties can thus be un- 
derstood from the coupling strength of the eigennodes. 
In harmonic systems, it is well known that 7^ is propor- 
tional to the square amplitude of the kth eigenmode at 
the chain end For fixed and free b.c, jk ~ k 2 /N 3 , 
7fe w l/N, respectively (for the sake of simplicity, we ne- 
glect the nonlinear dependence on fc, that is important 
only at large wavenumbers). Since the total heat flux is 
obtained by summing all contributions up to k — N, it 
follows that J is a quantity of order one independently 
of the b.c. 

The scenario becomes very different if some disorder is 
added, because the exponential localization of the eigen- 
modes makes their coupling with the heat baths basically 
negligible. In such conditions, only the first y/N modes 
are sufficiently extended to appreciably contribute to the 
heat flux. For free b.c, the total flux being the sum of 
the contribution is on the order of 1/ y^N which corre- 
sponds to a square root divergence of the condcutivity k. 
On the contrary, for fixed b.c, the sum of 7*, = k 2 /N 3 
yields J « 1/iV 3 / 2 , i.e. a vanishing conductivity. In other 
words, in the presence of disorder, a change of b.c. may 
turn a "heat superconductor" into a very good insulator!. 
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The question then arises of how the scenario modifies in 
the presence of nonlinearities. Past investigations of the 
FPU models clearly indicate that in the absence of inter- 
actions with thermal baths, the dynamics of a Fourier 
mode with small k is well described by the Langevin 
equation 



Qk+ulQk + ^ l Q k + Z nl =o 



(20) 



According to mode-coupling theory, j k ~ (k/N) 5 / 3 an 
estimate that is very close to the numerical values (see 
also Fig.Q. It thus appear natural to assume that in the 
presence of heat baths, the above equation modifies to 



Qk + to£Qk + (7* + 2jk)Qk + r + r + £ 



(21) 



although we anticipate that the effective value of the 
coupling 7^ with the heat bath cannot be the same as 
in the harmonic case. We start addressing the prob- 
lem of the effective temperature of the fcth mode. By 
repeating the same arguments that lead to predict a 
temperature T a in the harmonic case, we obtain now 
T k = [2 lk T a + J% l T£ l ]/(2 lk + Ik) wh ere we write T k 
to emphasize a possible dependence of the temperature 
on the wavenumber and where TJ} 1 is basically unknown, 
and it is expected to range between T~ and T + . For free 
b.c, 7fc ~ 1/-ZV, thus implying that there exists a criti- 
cal value k c ~ N 2 / 5 below which 7^ is negligible with 
respect to j k and viceversa above it. As a result, below 
k c , we expect T k = T a , while above, there should be a 
crossover towards the unknown value T}? 1 . Simulations 
performed in chains of different lengths with T + = 10 
and T~ = 8 are trivially in agreement with this scenario. 
Indeed, T k determined by evaluating the kinetic tempera- 
ture of the fcth mode is constantly equal to T a = 9. This 
implies that in spite of the nonlinearity in the profile, 
modal temperatures behave exactly as in the harmonic 
case (and also as when the Fourier law would apply). 
The same is true also for free boundary conditions when, 
being j k = k 2 /N 3 , there does not exist a value k c , since 
external dissipations should be always negligible. 

Let us now analyse the behaviour of modal fluxes. As 
the internal noise temperature adjusts itself to the arith- 
metic average of the heat-bath temperatures, the flux 
J^ 1 should vanish. This is confirmed by our simulations: 
the small values of J£ l that we have obtained can indeed 
be interpreted as statistical fluctuations. Much more in- 
triguing is the analysis of the true energy fluxes . With 
reference to free boundary conditions, the integral flux 
due to the modes with k < k c scales as I/A 3 / 5 and this 
would be consistent with the direct computations of heat 
conductivity, were the contribution of for k > k c neg- 
ligible. Unfortunately, if we assume that 7*, keeps scaling 
as 1/N (like in the harmonic case) we are led to conclude 
that all channels equally contribute to the heat flux with 
a consequent linear divergence of the conductivity with 
the system size! We already know that this is not the 
case, but the data reported in Fig. [S] give us some hints 
about possible explanations of this failure. In the figure 
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FIG. 8: Modal flux J„ (multiplied by the chain length N 
versus the scaled wavenumber k/N 2 ^ 5 for an FPU-/3 lattice 
with free boundary conditions. The solid, dashed and dotted 
lines refer to N = 256, 512 and 1024, respectively. Averages 
are taken over a time span of 10 8 units. 



we report the modal flux multiplied by N versus k/N 2 / 5 
for three different chain lengths. The reasonably good 
overlap confirms the hyphotesis that for k < k c w A 2 / 5 , 
the internal dissipation is negligible and the system be- 
haves like a harmonic chain. On the other hand, the 
decrease of J k suggests that above k c the effective cou- 
pling with the thermal baths is much smaller than 1/N. 
Determining the dependence of j k on k in this regime is, 
to our knowledge, an open problem. Solving this problem 
is all the way more crucial in the context of fixed bound- 
ary conditions where j k must differ everywhere from the 
harmonic chains. In fact, simulations publised elsewhere 
[f| indicate that, at variance with the harmonic case, al- 
though J k vanishes for k — > 0, the overall scaling be- 
haviour remains the same as in the free b.c. case. 

Finally, we wish to mention that the modal analysis 
can help to shed some light on a futher open problem: 
the determination of the temperature profile. We have 
seen that different models (such as harmonic and FPU 
chains) are characterized by energy equipartion even out 
of equilibrium. This is true even though they exhibit very 
different spatial profiles. It is quite natural to conjecture 
that this is caused by different phase correlations among 
the various modes, but working out a detailed explana- 
tion appears not to be a trivial task. 



VI. THE 2D FPU 

The presence of an energy threshold for the relaxation 
dynamics has been observed also in 2d lattice models 
of anharmonically coupled oscillators. For instance, in 
[40| the authors investigated a square-lattice of oscilla- 
tors coupled by a Lennard- Jones 6/12 potential 



V(r)=A (-) i2 -(-) 



(22) 



where r = \r\ is the modulus of the distance vector 
r between nearest-neighbor oscillators. At energy den- 
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sity e c « 0.3 they observed a crossover from a high- 
temperature regime of strong chaos to a low-temperature 
weakly chaotic regime. Specifically, above e c wave- 
packet excitations of low wave-number harmonic modes 
were found to relax very rapidly to a thermalized state, 
characterized by the equipartition of the energy among 
the Fourier modes. Below e c the weakly chaotic relax- 
ation dynamics exhibited a slowing-down of the energy 
equipartition, and the evolution appeared increasingly far 
from ergodic, while decreasing the energy density. In this 
regime, fluctuations of particles around their equilibrium 
positions are quite small, so that potential lll'L'l) can be 
very well approximated by its Taylor series expansion 

V(r) = ^ujr 2 + iar 3 + i/?r 4 + • • • (23) 

where the parameters to, a and (3 can be expressed in 
terms of A and a. This indicates that a very similar 
scenario is expected to hold also in the 2d FPU model, 
whose potential has the form reported on the r.h.s. of 
(EH- 

More recently, the FPU model has been investigated 
on a triangular 2d lattice |4l| • It has been observed that 
the slowing-down of energy equipartition for sufficienty 
small values of e is definitely less a dramatic effect than 
in the square lattice case. In particular, while increasing 
the size of the system there is evidence that in the tri- 
angular lattice a threshold value is better identified by 
the total energy E , rather than by the energy density e. 
This implies that the effects of weak chaos should van- 
ish in the thermodynamic limit much more rapidly than 
in the square-lattice case. Nonetheless, all of these re- 
sults can be interpreted consistently by considering that 
a typical feature of lattices of anharmonic oscillators is 
the appearance of long time scales for relaxation of exci- 
tations towards a thermalized state (equilibrium) as soon 
as the energy is sufficiently small. Accordingly, this fact 
is expected to have some important consequences on the 
heat transport also in 2d lattices. Actually, this prob- 
lem has been inv estig ated for the FPU and the Lennard- 
Jones models in |42j. Numerical simulations performed 
for lattice size N up to O(10 2 ) indicate that in the strong 
chaotic regime the thermal conductivity diverges loga- 
rithmically with the system size, n{N) ~ In TV. This is 
consistent with the theoretical prediction of the mode- 
coupling theory, which is obtained by estimating the av- 
erage value of the heat-flux time correlation function at 
thermal equilibrium jjj. Conversely, in the weak chaotic 
regime the thermal conductivity seems to diverge accord- 
ing to a power-law, k ~ N 71 . The exponent r\ increases 
while decreasing e below e c and it is expected to ap- 
proach unit for vanishing e, since in this limit the har- 
monic lattice is recovered. In analogy with what ob- 
served in Id, it seems reasonable to assume that also 
in 2d one should eventually recover the logarithmic di- 
vergence also in the weak chaotic regime for sufficiently 
large system size and sufficiently long integration time. 
On the other hand, the practical difficulty of performing 



more extended numerical simulations prevents the pos- 
sibility of verifying this conjecture. Moreover, it should 
be mentioned that there are also conflicting results indi- 
cating that the 2d scenario is still quite far from being 
fully understood. For instance, further numerical sim- 
ulations performed for the 2d FPU square-lattice were 
found to be compatible with a power-law divergence of 
the heat conductivity, with an exponent a « 0.22 |35|. 
In these simulations the authors used larger system sizes 
than those used in . Accordingly, their results should 
be considered more reliable for what the dependence of 
k on N is concerned, although noone among the avail- 
able theoretical arguments sup por ts such a prediction. 
It should be also noted that in [^a] the measurement of 
k(N) was obtained by non-equilibrium simulations with 
heat baths at relatively moderate temperatures with re- 
spect to those employed in ^3. Although the values of 
the temperature chosen in |35| are still in the range of 
strong chaos, they are not far from e c and one cannot 
exclude that the observed power-law behavior could be 
just an artifact of finite-size corrections already effective 
at temperatures too close to the so-called equipartition 
threshold. 



VII. CONCLUSIONS 

In spite of the many efforts made in the 50 years that 
separate us from the first numerical simulation performed 
by Fermi, Pasta and Ulam, one cannot yet onclude that 
heat conduction is fully understood. We discussed how 
the simple, apparently harmless, FPU model may serve 
to illustrate how both weakly chaotic dynamics and re- 
duced dimensionality affect the validity of macroscopic 
transport laws. 

Instances of open problems are the shape of the tem- 
perature profile which depends on boundary conditions 
and coexist with energy equipartion like in equilibrium, 
when temperature is constant across the system. The ef- 
fective coupling between thermal baths and (low- and 
high-frequency) Fourier modes is a related problem, 
whose solution might help in identifying those bound- 
ary conditions which can minimize the contact thermal 
resistance. A deeper understanding of these issues in 
toy models like FPU may be illuminating to describe en- 
ergy transport in "small systems" like single-molecules 
or nanostructured materials. 

The dependence of the time (and, correspondingly, the 
length) to reach the asymptotic scaling regime on the 
energy density at low temperatures represents another 
open problem. The results here reported indicate that 
the time needed for a nonlinear behaviour of the hydro- 
dynamic modes to set in increases as an inverse power of 
the energy density. The scaling behaviour is compatible 
with that exhibited by the times required for the chao- 
tization of generic trajectories, although it follows from 
the analysis of completely different dynamical processes. 
Accordingly, we are led to conjecture that in the thermo- 
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dynamic limit the evolution is controlled by a single time 
scale. Last but not least, the scenario in 2d FPU lattices 
is even more unclear, the nature of the leading behaviour 
being still questioned. 
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